Spreading with evaporation and condensation in one-component fluids 
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We investigate the dynamics of spreading of a small liquid droplet in gas in a one-component 
simple fluid, where the temperature is inhomogeneous around O.OTc and latent heat is released or 
generated at the interface upon evaporation or condensation (with Tc being the critical temperature). 
In the scheme of the dynamic van der Waals theory, the hydrodynamic equations containing the 
gradient stress are solved in the axisymmetric geometry. We assume that the substrate has a finite 
thickness and its temperature obeys the thermal diffusion equation. A precursor film then spreads 
ahead of the bulk droplet itself in the complete wetting condition. Cooling the substrate enhances 
condensation of gas onto the advancing film, which mostly takes place near the film edge and can be 
the dominant mechanism of the film growth in a late stage. The generated latent heat produces a 
temperature peak or a hot spot in the gas region near the film edge. On the other hand, heating the 
substrate induces evaporation all over the interface. For weak heating, a steady-state circular thin 
film can be formed on the substrate. For stronger heating, evaporation dominates over condensation, 
leading to eventual disappearance of the liquid region. 

PACS numbers: 68.03. Fg, 68.08.Bc, 44.35.-l-c, 64.70.F- 
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I. INTRODUCTION 

Extensive efforts have been made on the static and 
dynamic properties of wetting transitions for various flu- 
ids and substrates both theoretically and experimentally 
[l|. In particular, spreading of a liquid has been stud- 
ied by many groups [l|-l5| , since it is of great importance 
in a number of practical situations such as lubrication, 
adhesion, and painting. Hydrodynamic theories were de- 
veloped for spreading of an involatile liquid droplet in 
gas in an early stage of the theoretical research [ll, S] • A 
unique feature revealed by experiments [2, ^ |7| is that a 
thin precursor fllm is formed ahead of the liquid droplet 
itself in the complete wetting condition. Hardy flrst re- 
ported its formation ascribing its origin to condensation 
at the film edge Jal, but it has been observed also for 
involatile fluids [g, 0] ■ To understand nanometer-scale 
spreading processes, a number of microscopic simulations 
have been performed mainly for fluids composed of chain- 
like molecules fsl-flSi. 



However, understanding of the wetting dynamics of 
volatile liquids is still inadequate. We mention some ex- 
amples where evaporation and condensation come into 
play. In their molecular dynamic simulation [13], Koplik 
et al. observed evaporation of a droplet and a decrease of 
the contact angle upon heating a substrate in the partial 
wetting condition. In their experiment [16|, Guena et al. 
observed that a weakly volatile droplet spread as an in- 
volatile droplet in an initial stage but disappeared after 
a long time due to evaporation in the complete wetting 
condition. In a near-critical one-component fluid [iTj, 
Hcgscth et al. observed that a bubble was attracted to a 
heated wall even when it was completely wetted by liq- 
uid in equilibrium (at zero heat flux), where the apparent 
contact angle of a bubble increased with the heat flux. 

In addition to spreading on a heated or cooled sub- 
strate, there are a variety of situations such as droplet 



evaporation |18l - l22j . boihng on a heated substrate l23l- 
|25| , and motion of a bubble suspended in liquid [26|, |23l , 
where latent heat generated or released at the interface 
drastically influences the hydrodynamic processes. In 
particular, a large temperature gradient and a large heat 
flux should be produced around the edge of a liquid film 
or the contact line of a droplet or bubble on a substrate 
[23 . |24| . The temperature and velocity profiles should be 
highly singular in these narrow regions. Here an experi- 
ment by Hohmann and Stephan [2g| is noteworthy. They 
observed a sharp drop in the substrate temperature near 
the contact line of a growing bubble in boiling. Further- 
more, we should stress relevance of the Marangoni flow 
in multi-component fluids in two-phase hydrodynamics 
0, [2J, 123] , where temperature and concentration vari- 
ations cause a surface tension gradient and a balance of 
the tangential stress induces a flow on the droplet scale. 

In hydrodynamic theories, the gas-liquid transition has 
been included with the aid of a phenomenological input 
of the evaporation rate on the interface J . Some authors 
[l8l - [20 ] assumed the form J{r, t) = Ja/\/re{t)^ — r^ for a 
thin circular droplet as a function of the distance r from 
the droplet center, where re(i) is the film radius and Jg 
is a constant. In the framework of the lubrication the- 
ory, Anderson and Dabis |29| examined spreading of a 
thin volatile droplet on a heated substrate by assuming 
the form J = {Tj ~ Tc^)/K* , where T/ is the interface 
temperature, Tex is the saturation (coexistence) temper- 
ature, and K* is a kinetic coefficient. In these papers, 
the dynamical processes in the gas have been neglected. 

Various mesoscopic (coarse-grained) simulation meth- 
ods have also been used to investigate two-fluid hydro- 
dynamics, where the interface has a finite thickness. We 
mention phase field models of fiuids (mostly treating in- 
compressible binary mixtures) j3C|-i45j] , where the gradi- 
ent stress is included in the hydrodynamic equations (see 
a review in Ref.[32(). In particular, some authors numer- 



ically studied liquid-liquid phase separation in heat flow 
3lL aa . l40l . l44| , but these authors treated symmetric bi- 
nary mixtures without latent heat. Recently, one of the 
present authors developed a phase field model for com- 
pressible fluids with inhomogeneous temperature, which 
is called the dynamic van der Waals model [4l|, 1421 . In its 
framework, we may describe gas-liquid transitions and 
convective latent heat transport without assuming any 
evaporation formula. In one of its applications [22| . it 
was used to investigate evaporation of an axisymmetric 
droplet on a heated substrate in a one-component system. 
Our finding there is that evaporation occurs mostly near 
the contact line. We also mention the lattice Boltzmann 
method to simulate the continuum equations, where the 
molecular velocity takes discrete values [ST'-HlO', 45] . How- 
ever, this method has not yet been fully developed to 
describe evaporation and condensation. 

In this paper, we will simulate spreading using the dy- 
namic van der Waals model |,41, 42]. We will treat a one- 
component fluid in a temperature range around 0.9Tc, 
where the gas and liquid densities are not much sepa- 
rated. Namely, we will approach the problem relatively 
close to the critical point. Then the mean free path in the 
gas is not long, so that the temerature may be treated 
to be continuous across an interface in nonequilibrium. 
When the gas is dilute, the phase field aproach becomes 
more difficult to treat gas flow produced by evaporation 
and condenssation. It is known that the temperature 
near an interface changes sharply in the gas over the 
mean free path during evaporation [46J. 

The organization of this paper is as follows. In Sec. II, 
we will present the dynamic equations with appropriate 
boundary conditions. In Sec. Ill, the simulation method 
will be explained. In Sec. IV, numerical results of spread- 
ing will be given for cooling and heating the substrate. 



II. DYNAMIC VAN DER WAALS THEORY 

When we discuss phase transitions with inhomoge- 
neous temperature, the free energy functional is not well 
defined. In such cases, we should start with an entropy 
functional including a gradient contribution, which is de- 
termined by the number density n = n(r,t) and the in- 
ternal energy density e = e(r, t) in one-component fluids. 
Here we present minimal forms of the entropy functional 
and the dynamic equations needed for our simulation. 



depend on n [42] , but it will be assumed to be a positive 
constant independent of n. The gradient entropy is nega- 
tive and is particularly important in the interface region. 
The entropy functional is the space integral Sb = J drS 
in the bulk region. As a function of n and e, the temper- 
ature T is determined from 
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(2.2) 



The generalized chemical potential /t including the gra- 
dient part is of the form. 



\ on 



fi - TCV^n, 



(2.3) 



where /i = —T[d{ns)/dn\e is the usual chemical poten- 
tial per particle. In equilibrium T and fi are homoge- 
neous constants. In this paper, we introduce the gradi- 
ent entropy as in Eq.(2.1), neglecting the gradient energy 
[4l|,|42l . Then the total internal energy in the bulk is sim- 
ply the integral / dre. 

In the van der Waals theory [43], fluids are character- 
ized by the molecular volume vq and the pair-interaction 
energy e. As a function of n and e, s is written as 



== ks ln[(e/7i 4- evo72)^/'^{l/von ~ 1)] -|- 



so, 



(2.4) 



where SQJks ~ ln]i'o(TO/37rfi,^)'^/^] -I- 5/2 with to being 
the molecular mass. We define T as in Eq.(2.2) to obtain 
the well-known expression for the internal energy e — 
inksT /2 — evQ'n? and the pressure 

p — n^ + TS — e = nkBT/{l — Vqu) — evQri^ . (2.5) 

The critical density, temperature, and pressure read 

ne = l/3i;o, Te = 8e/27fcs, p^ = e/27va, (2.6) 

respectively. Macroscopic gas-liquid coexistence with a 
planar interface is realized for T < T,. and at the satu- 
rated vapor pressure p = PcxiT). With introduction of 
the gradient entropy, there arises a length (. defined by 



i={Cl2kBVof'\ 



in addition to the molecular diameter 



,1/3 



(2.7) 



From 



A. Entropy formalism 



Eq.(2.3) the correlation length £^ is defined by ,J ^ 
{dy,/dn)T/TC^ so ^ is proportional to £ as 



We introduce a local entropy density S = S(r,t) con- 



sisting of regular and gradient terms as [4l|, |42| 



S = ns{n,e)~-C\'^n\^. 



(2.1) 



Here s = s{r^t) is the entropy per particle depending 
on n and e. The coefficient C of the gradient term can 



^/£ = ni2vokBTKT)^^^, 



(2., 



where Kt = {dn / dp)T / n is the isothermal compressibil- 
ity. The interface thickness is of order £^ in two-phase 
coexistence. The ratio ^/wq should be of order unity 
for real simple fiuids. However, we may treat £ as an 
arbitrary parameter in our phase field scheme. 



B. Hydrodynamic equations 



C. Boundary conditions 



We set up the hydrodynamic equations from the princi- 
ple of positive entropy production in nonequihbrium [48J . 
The mass density p = mn obeys the continuity equation, 






(2.9) 



where v is the velocity field assumed to vanish on all the 
boundaries. In the presence of an externally applied po- 
tential field U{r) (per unit mass), we write the equation 
for the momentum density pv as 



di 



pv = —V • {pvv + n — cr) — pWU. 



(2.10) 



In our previous work[42j we set U — gz for a gravitational 
field with g being the gravity acceleration. We note that 
U may also represent the van der Waals interaction be- 
tween the fluid particles and the solid depending the dis- 
tance from the wall [l| . The stress tensor is divided into 

three parts. The pvv is the inertial part. The 11 — {11^} 
is the reversible part including the gradient stress tensor. 



n,; 



p-CT(nV2n+-|V7ip) 



(2.11) 



where p is the van der Waals pressure in Eq.(2.5). Here- 
after Vi = d/dxi with Xi representing x, y, or z. The 
a = {aij} is the viscous stress tensor expressed as 

cTrj = ry(V,Wj- + Vji-O + (C - 277/3)(V • v)6^j, (2.12) 

in terms of the shear viscosity rj and the bulk viscosity 
C- Including the kinetic energy density and the potential 
energy, we define the (total) energy density by ex = e -I- 
pv'^/2 + pU. It is a conserved quantity governed by [43] 



d_ 
di 



ex = —V • 



eTV + {Il-a)-v- AVT 



(2.13) 



where A is the thermal conductivity. With these hydro- 
dynamic equations including the gradient contributions, 
the entropy density S in Eq.(l) obeys 






Sv-Cn{\7.v)\/n--\/T 



(e 



T 



, (2.14) 



where the right hand side is the nonnegative-definite en- 
tropy production rate with 



ei, = 2J TijVjVi, 60 = A(VT)^/r. 



(2.15) 



In passing, the constant Sq in Eq.(2.4) may be omitted 
in Eq.(2.14) owing to the continuity equation (2.9). 



We assume the no-slip boundary condition, 

v = 0, (2.16) 

on all the boundaries for simplicity. However, a number 
of molecular dynamic simulations have shown that a slip 
of the tangential fluid velocity becomes significant around 
a moving contact line [50, [Sll . 

We assume the surface entropy density (Ts('t-s) and the 
surface energy density £5(^5) depending on the fluid den- 
sity at the surface, written as n^. The total entropy in- 
cluding the surface contribution is of the form. 



drS 



daa^ 



(2.17) 



where / da is the surface integral over the boundaries. 
The total fluid energy is given by 

1 



dr{e + -pv^ + pU) + / dae 



(2.18) 



We assume that there is no strong adsorption of the 
fluid particles onto the boundary walls. The fluid den- 
sity is continuously connected from the bulk to the 
boundary surfaces; for example, we have ns{x,y,t) — 
limz_>_|_o ".(r, i) at z = 0. Then the total particle number 
of the fluid in the cell is the bulk integral N = j dm. 

We assume that the temperatures in the fluid and in 
the solid are continuously connected at the surfaces. The 
temperature on the substrate is then well-defined and we 
may introduce the surface Helmholtz free energy density, 

f, = es-Tas. (2.19) 

As the surface boundary condition, we require 



TCub ■ Vn 



(2.20) 



where Oh is the outward surface normal unit vector. 
This boundary condition has been obtained in equilib- 
rium with homogeneous T by minimization of the total 
Helmholtz (Ginzburg-Landau) free energy, 

i^tot = / dr{e -TS)+ f daf,. (2.21) 

We assume this boundary condition in Eq.(2.20) even in 
nonequihbrium. Then use of Eq.(2.14) yields 48] 

Ob ■ AVT + e 



^5 - 

-TT'-'tot — 

dt 



dr- 



ee 



T 



da- 



T 



(2.22) 



where Cs — dcs/dt = {des / dns){dns / dt) . The first term 
in the right hand side is the bulk entropy production rate, 
while the second term is the the surface integral of the 
heat flux from the solid divided by T or the entropy input 
from the solid to the fluid. 

In this paper, we present simulation results with [7 = 
for simplicity. In our previous work [42| a large gravity 
field was assumed in boiling. In future we should investi- 
gate the effect of the long-range van der Waals interaction 
in the wetting dynamics. 



III. SIMULATION METHOD 

In our phase field simulation, we integrated the con- 
tinuity equation (2.9), the momentum equation (2.10), 
and the entropy equation (2.14), not using the energy 
equation (2.13), as in our previous simulation [23. With 
this method, if there is no applied heat flow, temperature 
and velocity gradients tend to vanish at long times in the 
whole space including the interface region. This numeri- 
cal stability is achieved because the heat production rate 
iv + eg > appears explicitly in the entropy equation, 
so that dStot/dt > in Eq.(2.22) without applied heat 
flow. We can thus successfully describe the temperature 
and velocity near the film edge (those around the contact 
line of an evaporating droplet in Ref.|2a|). 

It is worth noting that many authors have encountered 
a parasitic flow around a curved interface in numerically 
solving the hydrodynamic equations in two-phase states 
[45I |52| . It remains nonvanishing even when the system 
should tend to equilibrium without applied heat flow. It 
is an artificial flow, since its magnitude depends on the 
discretization method. 



The transport coefRcients are proportional to cr^/^. In 
this paper we set a = 0.06, for which sound waves are 
well-defined as oscillatory modes for wavelengths longer 
than^ (see Fig. 5) [13]. 

The temperature at the top z = H is fixed at Th, 
while the side wall at r = L is thermally insulating or 
I'h • VT = dT/dr = at r = L. The boundary condition 
of the density n on the substrate z = is given by 



voi 



dn 

dz 



-$i. 



(3.4) 



where $1 arises from the short-range interaction between 
the fluid and the solid wall [l|? ]. We treat $1 as a 
parameter independent of T. From Eq.(2.19) this can be 
the case where es = and as{ns) = {^^iC /v^tjUs- For 
example, at T = 0.875Tj,, the contact angle is zero at 
^i ^ 0.060 and the wall is completely wetted by liquid 
for larger $1. Furthermore, we set dn/dz = on the top 
plate at z ^ H and dn/dr = on the side wall at r = L. 



Solid substrate 



A. Fluid in a cylindrical cell 

We suppose a cylindrical cell. Our model fluid is in 
the region < z < H and < r = (x^ + y^)^/^ < ^^ 
where H = 300Ax and L = 400Aa; with Ax being the 
simulation mesh length. The velocity field v vanishes on 
all the boundaries. In this axisymmetric geometry, all the 
variables are assumed to depend only on z, r and t. The 
integration of the dynamic equations is on a 200 x 400 
lattice in the fluid region. We set Ax — £/2, where £ is 
defined in Eq.(2.7). We will measure space in units of i. 
Then H = 150 and L = 200 in units of £. 

The transport coefficients are proportional to n as 



rj = C, = VQmn, A = ksi^QTi. 



(3.1) 



These coefficients are larger in liquid than in gas by the 
density ratio ni/ng{'^ 5 in our simulation). The kine- 
matic viscosity vq = rj/mn is a constant. We will mea- 
sure time in units of the viscous relaxation time. 



To = f/i^o = C/2kBVaiyo, 



(3.2) 



on the scale of £. We will measure velocities in units of 
£/to = Vol £. The time mesh size of our simulation is 
Ai = O.OItq. Away from the criticality, the thermal dif- 
fusivity Dt = ^/Cp is of order vq and and the Prandtl 
number Pr ~ vq/Dt is of order unity, so tq is also the 
thermal relaxation time on the scale of £. Here the iso- 
baric specific heat Cp per unit volume is of order n far 
from the criticality, while it grows in its vicinity. With 
Eq.(3.1), there arises a dimensionless number given by 



In our previous work, we assumed a constant temper- 
ature at the bottom plate z = [23 . [4ll . [43 |. In this 
paper, we suppose the presence of a sold wall in the re- 
gion — iJ^, < z < and < r = {x^ + y'^Y^^ < L, where 
its thickness is i/^, = lOOAx = 50-^ = H/3. The temper- 
ature in the solid obeys the thermal diffusion equation. 



^t«~:Tr — Atu V 1 , 
at 



(3.5) 



where C^ is the heat capacity (per unit volume) and A^ 
is the thermal conductivity of the solid. The temperature 
T{r, z, t) is continuous across the substrate z = 0. In our 
simulation, the thermal diffusivity in the solid is given 
by Dw — Xw/Cw — 400i'o, while the thermal diffusivity 
of the fluid Dt is of order vq away from the critical- 
ity. Thus the thermal relaxation time in the substrate 
is H^/Dyj = 25ro, which is shorter than typical spread- 
ing times to follow. Because D^ 3> Dt, we integrated 
Eq.(3.5) using the implicit Crank-Nicolson method on a 
100 X 400 lattice. 

In this paper, the temperature T at the substrate bot- 
tom z — —Hw is held fixed at a constant Tw That is, 
for any r, we assume 



T{r, -7J„) = T^ 



(3.6) 



Heating (cooling) of the fluid occurs when r„, is higher 
(lower) than the initial fluid temperature Tq. There is no 
heat flux through the side wall, so dT/dr — Q at r — L 
as in the fluid region. From the energy conservation at 
the boundary, the heat flux on the substrate surface is 
continuous as 



mvl/e£^ = m£^/eTl. 



(3.3) 



(A.,r') 



{XT') 



=4-0, 



(3.7) 



where T' = dT/dz. This holds if there is no appreciable 
variation of the surface energy density Es . We define the 
parameter, 



A = A/(rwoA„,) = ksva/v^K 



(3.8) 



Then {T')z=_q — AnsUo(r')z=+o on the substrate. In 
this paper, A is set equal to 0.002 or 0.2. We found that 
the boundary temperature at z = is nearly isothermal 
'AiT — Tw for A — 0.002 but considerably inhomogeneous 
around the edge for A = 0.2. 



C. Preparation of the initial state and weak 
adsorption preexisting before spreading 

To prepare the initial state, we first placed a semi- 
spheric liquid droplet with radius R — A0£ on the sub- 
strate z = with gas surrounding it. Here we set $i = 
to suppress adsorption of the fluid to the solid. The 
temperature and pressure were T — Tq ~ 0.875Tc and 
P = Pcx{To) = 0.573pc on the coexistence line in the fluid. 
The liquid and gas densities were those on the coexis- 
tence curve, n° — 0.579W{^^ in liquid and n^ — 0.123w^^ 
in gas. The entropy difference between the two phases 
is 2.1kB per particle. The total particle number is 
N = 27r(n^ - n°g)R^/3 + nn^gL^H = 1.61 x lOS^s/^;^. 
The particle number in the droplet is about 5% oi N. 

Next, we waited for an equilibration time of 10^ with 
$1 — 0. The contact angle was kept at 7r/2 and 
Ob ■ Vn = on all the boundary surfaces. However, 
the liquid and gas pressures were slightly changed to 
0.608pc and 0.575pc, respectively. The pressure differ- 
ence Ap = 0.033pc is equal to 2^/R from the Laplace 
law. In accord with this, the surface tension 7 at T = Tq 
is given by 7 = 0.66ipc in our model. As a result, the 
liquid density was increased to O.SSSu,^^ and the droplet 
radius was decreased to 38€. After this equilibration we 
hereafter set t = as the origin of the time axis. 

At t = 0, we changed the wetting parameter $1 in the 
boundary condition (3.4) from to 0.0610 to realize the 
complete wetting condition. Before appreciable spread- 
ing, weak adsorption of the fluid has been induced on 
the substrate in a short time of order unity (in units of 
To). For small $1 and away from the contact line, this 
preexisting density deviation, written as 6n{z), is of the 
exponential form. 



Sn{z) = (^$1/1-0^)6 



-^/« 



(3.9) 



in terms of the correlation length ^. Note that homo- 
geneity of fi in Eq.(2.3) yields (f^^ — d'^/dz^)Sn = in 
the linear order, leading to Eq.(3.9) under Eq.(3.4). The 
z integration of Sn(z) is the excess adsorption. 



ad 



e'*iM^- 



(3.10) 



In the gas at T == 0.875Tc, Eq.(2.8) gives $ = 1.68£, lead- 
ing to Fad — 0.24^/wo. We shall see that this adsorption 
is one order of magnitude smaller than that due to a 
precursor fllm (^ 2.51/ vq in Fig. 6 below). 



IV. SPREADING ON A COOLED SUBSTRATE 

We present numerical results of droplet spreading on 
a cooler substrate. At i = the bottom temperature T^ 
at z = -Hy, was lowered from Tq = 0.875rc to 0.870rc 
except for two curves in Fig. 2 (for which T^, = Tq even for 
t > 0)). The top temperature at z = H was kept at Tq in 
all the cases. Subsequently, we observed spreading with 
an increase of the liquid fraction due to condensation. 
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FIG. 1: (Color online) Shapes of an axisymmetric droplet 
spreading on a cooler substrate with T^ — 0.870Tc at various 
times for A = 0.002 (top) and 0.2 (middle) in the r-z plane. 
The system temperature was initially To = 0.875Tc at i = 0. 
The boundary position between the main body of the droplet 
and the precursor film is fixed at r = rth = 52.5^ for both 
A. The edge position re{t) of the film increases with time as 
illustrated in the bottom plates. 
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FIG. 2: (Color online) Time evolutions of the edge position 
re(i) divided by L (left) and the particle number in the droplet 
Ne{t) divided by Ne{0) (right) for A = 0.2. Two curves cor- 
respond to r„ = 0.870Tc < To (red) and T^ = 0.875rc = To 
(green). The interface curve is determined by Eq.(4.2). The 
film edge reaches the side wall at i ~ 10*. Condensation 
occurs faster in the cooled case than in the non-cooled case. 
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FIG. 3: (Color online) Heat flux on the substrate Q]^[r,t) 
as a function of r in units of el/vQTo at various times for 
r„ = O.STOTc < To with A = 0.002 and 0.2. A negative 
peak at the film edge indicates absorption of latent heat from 
the fiuid to the solid. At long times this absorption becomes 
vifeaker and there appears a heat flow from the solid to the 
fluid for r < rth. 



A. Evolution on long and short time scales 

In Fig.l, the droplet spreads over the substrate in 
the complete wetting condition for A — 0.002 and 0.2. 
The liquid region is divided into the droplet body in 
the region r < r^h and the precursor film in the re- 
gion rth < r < feit). In our simulation, rth is equal 
to 52.5^ = 0.26L independently of time, while re(i) in- 
creased in time. The film thickness £f was only weakly 
dependent on time being about 5£ for both A (see the 
film profiles in Fig. 6 below). However, for slightly deeper 
cooling (say, for T„ = 0.868Tc) or for slightly larger $i 
(say, for $i = 0.065), a new liquid region (a ring here) 
appeared on the substrate ahead of the precursor film. 

In Fig. 2, we show re(t) and the particle number in 
the liquid region N£{t) vs i for A = 0.2 in the cooled 
case with T^ = O.STOTc and the non-cooled case with 
r,„ = To = 0.875Tc. We calculate Neit) from 



N,{t) - 27r 



= («) 



drr 



Zint (»•,*) 



dz n{r,t), (4.1) 



where the interface height is at z = 2int(?', t) in the range 
< r < re(t). It starts from the initial number iV^(O) = 
0.67x 10^^^/wo and becomes a few times larger at t ^ 10**. 
Here condensation takes place even for the non-cooled 
case with T^ = Tq. In these two cases, the latent heat due 
to condensation is mostly absorbed by the solid reservoir. 
In calculating Ni(t) we determine the film height zij^tir, t) 
from the relation, 



n(r,2int,i) = (n^-Kn°)/2, 



(4.2) 



where n^ = 0.579uq and n" 



0.1 23^0^^ are the densities 
on the coexistence curve at T = O.STSTc. In our case, the 
film is so thin and there is no unique definition of Zint- 
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FIG. 4: (Color online) Temperature T around an advancing 
fllm edge at t = 1000, where T™ = 0.870Tc and A = 0.2. In 
the top, the color represents the temperature according to the 
color map, and the velocity fleld is shown by arrows with its 
maximum being 1.4 x 10~^i'/ro as indicated below the plate. 
In the bottom, the substrate temperature at z = is plotted, 
which is maximum at the edge position due to a flnite thermal 
conductivity of the solid. 
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FIG. 5: Pressure (left) and temperature (right) vs t at the 
point {z/ H, r/L) = (0.48, 0.5) far from the substrate in gas, 
where T^ — 0.870rc and A — 0.2. Their short time behavior 
(i < 200) is due to propagation of a low-pressure sound pulse 
and is adiabatic (inset), while their long time behavior is due 
to gradual condensation. 



In Fig. 3, we display the heat flux on the substrate 
Qh{r,t) for the same runs. From Eq.(3.7) it is defined 
in terms of the temperature gradient T' — dT/dz as 

Qb(r, t) - -(A^T').=-o = -(AT').=+o. (4.3) 

Negative peaks indicate absorption of latent heat from 
the fiuid to the substrate around the film edge. However, 
at long times {t = 5000 in the figure) heat is from the 
solid to the fluid in the region of the droplet body r < rth. 



The amplitude of Qij(r,t) around the peak is larger for 
A = 0.002 than for A — 0.2, obviously because heat is 
more quickly transported for smaller A or for larger A^;. 
Also Qh{r, t) is sensitive to Tq — T^. For example, in the 
non-cooled case T^ = Tq, the minima of Qhi^, t) became 
about half of those in Fig. 3 (not shown here). In our 
previous simulation [22], a positive peak of Qh{r,t) was 
found at the contact line of an evaporating droplet. 

In the upper panel of Fig. 4, we show the temperature 
near the edge a.t t — 1000, where A = 0.2 and T^ = 
0.870Tc < Tq. It exhibits a hot spot in the gas side 
produced by latent heat. In this run, the peak height of 
the hot spot Tp depended on t as 10'^ {Tp-Tyj)/Tc = 1.0, 
0.7,0.5, and 0.4 for lO^^i = 1, 2, 3, and 4. The maximum 
of the gas velocity Vg is 0.014 around the hot spot, while 
the edge speed is a few times faster as dre/dt ~ 0.04. 
The corresponding Reynolds number Vg£f/i>o in the gas is 
very small (^ 0.07 here). In the non-cooling case T^^ = Tq 
the peak height was reduced to Tp — Tq — 0.007rc and 
Vg to 0.008 at t = 10^. In the lower panel of Fig. 4, the 
substrate temperature at z = is maximum at the film 
edge. Such a temperature variation in the solid should 
be measurable [25|. 

In Fig. 5, we display the time evolution of the pres- 
sure and the temperature at the position {z,r) — 
{OASH, 0.5L) in the gas region far from the substrate 
in the case T^ = 0.870Tc and A = 0.2. In the inset, their 
initial deviations originate from a lower-pressure sound 
pulse emitted from the adsorption layer in Eq. (3.9). This 
acoustic process is an example of the piston effect |53l . l54| . 
In this case the thermal diffusion layer due to cooling of 
the substrate gives rise to a smaller effect. The emitted 
pulse traverses the cell on the acoustic time H/cg ^ 50 
and is reflected at the top plate, where Cg ^ 4 is the sound 
velocity in the gas. The deep minimum of T below T^ 
and that of p at t '^ 25 are due to its first passage. Here 
the adiabatic relation ST = {dT / dp) sSp is well satisfied 
for the deviations 6T = T ~ Tq and Sp = p — pq. The 
adiabatic coefficient [dT /dp)s is equal to llTc/pc in the 
gas and is larger than that in the liquid by one order of 
magnitude. On long time scales. Fig. 5 shows that the 
pressure gradually decreases with progress of condensa- 
tion, while the temperature increases for 200 ^ t < 1500, 
slowly decreases for 1500 < t ^ 3000, and again increases 
for longer t. The gas temperature in the middle region is 
shghtly higher than T„ by 0.002Tc at t = 9000. We note 
that the gas temperature is influenced by a gas flow from 
the droplet and behaves in a complicated manner. 



B. Profiles of density, temperature, and pressure 
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FIG. 6: Density (top), temperature (middle), and normal 
pressure (bottom) as functions of the distance z from the sub- 
strate at r/L — 0.125 (left) and at r/L — 0.55 (right), where 
t = 2000, r„ = O.STOTc, and A = 0.2. The former path passes 
through the droplet body, while the latter through the film 
edge. The black dot • on each curve indicates the interface 
position determined by Eq.(4.2). 
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FIG. 7: (Color online) Temperature T around a film edge in 
the fluid (red) and in the solid (black) in the r-z plane at 
t = 2000, where T™ = 0.870Tc and A = 0.2. See Fig.4 for the 
hot spot at f = 1000 in color in the same run. 



In Fig. 6, we show the profiles of the density n, the 
temperature T, and the the stress component p along 
the density gradient at i = 2000 for A — 0.2 and T^j — 
O.STOTc. We define p as 



P 



Z^ J>»i>jn,j 



p-CT{nV^n-\Vnf/2), (4.4) 



where 11^ is the reversible stress tensor in Eq.(2.11), p is 
the van der Waals pressure, and £> = {vi} — {Vin/|Vn|} 
is the unit vector along the density gradient Vn. Thus 
p is called the normal pressure. Obviously, p = p in the 
bulk region. In equilibrium, p is equal to the satura- 
tion pressure pcx(T) for a planar interface [42l |. while it 



changes by the Laplace pressure difference 27/i? along i> 
across an interface with mean curvature 1/ R. In nonequi- 
librium, we find that inhomogeneity of p around an in- 
terface is much weaker than that of p itself. The left 
panels for r = 0.125L < rth in Fig. 6 indicate weak ad- 
sorption near the wall in Eq.(3.9), a well-defined interface 
at z ^ 20, and a negative temperature gradient within 
the droplet body. For this r, a heat flow is from the solid 
to the fluid. In the right panels for r = 0.55L > rth in 
Fig. 6, n decreases from a liquid density near the wall to 
a gas density without a region of a flat density and T 
exhibits a peak at the hot spot. Furthermore, Fig. 7 gives 
a bird view of the temperature near the edge from the 
same run, which corresponds to the middle right panel in 
Fig. 6. Here the temperature inhomogeneity in the solid 
can also be seen. 

It is of interest how the normal pressure and the tem- 
perature {p,T) at the interface is close to the coexistence 
line {pcx{T),T) in the p-T phase diagram. We define 



h^ 



T-Tq 

Tr. 



dT \ p-po 



(4.5) 



where the derivative {dT/dp)cx along the coexistence line 
is equal to 0.38Tc/pc at T = 0.875Tc. The upper panel 
of Fig. 8 displays h around the film at i = 1000, while the 
lower panel of Fig. 8 gives h along the surface z = Zi^t at 
four times for A = 0.2 and 0.002. This quantity repre- 
sents the distance from the coexistence line p = Pcx{T). 
In the bulk region, h < m stable liquid and metastable 
gas, while ft. > in stable gas and metastable liquid. We 
can see that h nearly vanishes in the droplet body r < rth 
and increases in the film rth < r < ^eit), but h remains 
less than 10~^ even at the edge. Note that the Laplace 
pressure contribution to h is {dT / dp)cx'2^ /TcR, which is 
of order 0.01 in the droplet body r < rth at i = 1000. 



C. Condensation rate and gas velocity 

In our previous simulation|22|. evaporation of a thick 
liquid droplet mostly takes place in the vicinity of the 
contact line in the partial wetting condition. We here 
examine the space dependence of the condensation rate 
of a thin fim in the complete wetting condition. 

We introduce the number flux J{r, t) from gas to liquid 
along v — |Vri|^^Vn through the interface. 



J{r,t) = n{v -Vint) 



(4.6) 



where Vint is the interface velocity. If J is regarded as a 
function of the coordinate along the normal direction O, 
it is continuous through the interface from the number 
conservation, while n and v ■ i> change discontinuously. 
Thus we may well determine J on the interface. If it 
is positive, it represents the local condensation rate per 
unit area. In Fig. 9, we plot J{r, t) vs r / L in the region 
< r < re(i) at three times for A = 0.002 and 0.2 in 
the cooled case T^ = 0.870Tc. We recognize that J{r, t) 




0.1 0.2 0,3 0.4 0.5 

-0,006 0012 



A=0.2 



(a) 10 h 



ai (a) 




t=1000 
t=2000 _ 
t=3000 
t=4000 



0.4 

r/L 



FIG. 8: (Color online) Top: Distance from the coexistence 
line h in Eq.(4.5) a.t t = 1000 in color, which is negative in 
the droplet body and is positive in the film and in the gas. 
Here T„ = O.STOTc and A = 0.2. Bottom: (a) h and (b) 
(T — To)/Tc along the interface at four times. For r < rth, \h\ 
is smaller than (To — T)/Tc. For r > rth, the distance from 
the coexistence line increases. 



steeply increases in the precursor film and is maximum 
at the edge. Moreover, it becomes negative in the body 
part r < rth at t = 3000, where evaporation occurs. 
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FIG. 9: (Color online) Flux J(r, t) on the interface in units 
of I/vqtq vs r/L in the region < r < re{t) at 10~''i = 1, 
2, and 3 for T^ = O.STOTc in the two cases of A = 0.002 and 
0.2. A precursor film is on the left of the arrow (see Fig.l). In 
its positive region it is the condensation rate. In its negative 
region evaporation takes place. 



The total condensation rate Wtot{t) is the surface in- 
tegral of J{r^ t) on all the surface. The surface area in 
the range [r,r + dr\ is da — 2TTdrr/ sind, where 9 is the 
angle between u and the r axis. Thus, 



Wtotit) = 2TT dr rJ{r,t) /sin ( 



(4.7) 
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FIG. 10: (Color online) Total condensation rate Wtotit) 
(green), condensation rate onto the film Wfiim(i) (red), and 
flow from the droplet body to the film Jflow(i) (blue) in 
units of £^/voTo as functions of time. The time range is 
[3 X 10^7 X 10^] for A = 0.002 (left) and [3.3 x 10^ 7 x 10^] 
for A = 0.2 (right). 



The particle number in the liquid region Ni{t) in Eq.(4.1) 
increases in time as 



-Ne{t) = Wtotit). 
at 



(4.8) 



We also define the condensation rate in the film region, 

Wiiim{t) = 2Tr drrj{r,t)/sm0, (4.9) 

where sin = 1. In this integral the vicinity of the edge 
gives rise to a main contribution. In fact, the contribu- 
tion from the region re — 16£ < r < re is about 50% of the 
total contribution from the region rth < r < r^. There- 
fore, in terms of the gas velocity Vg and the gas density 
Ug around the edge, we estimate Wfiim{t) as 



Wfiunit) -- 2TrrengVg£c, 



(4.10) 



where £c is the width of the condensation area estimated 
to be about 30^. 

The flux from the droplet body to the film is given by 

Jflow (i) = 27rrth / dz n{rth, z,t)vr{rth, z,t), (4.11) 
Jo 

where Vr{r,z,t) = v^x/r + Vyy/r is the velocity in the 
plane within the film. This lateral flux is defined at r = 
rth- More generally, we may introduce the flux 

Jf(r, t) = 27rr / dz n{r, z,t)vr{r, z,t), (4.12) 

Jo 

for r > rth with zth being the film thickness £f. Then 
Jfiowit) = J{{rth,t). In the presence of condensation onto 
the film, J{{r,t) increases with increasing r. Its max- 
imum Jf{re{t),t) is estimated as 2TrrefiiV££f, where Hi 
is the average density in the film and vg is the average 
fiuid velocity in the film. At i = 1000 (or 4000), the ratio 
Jf(re(i),t)/Jf(rth,i) was equal to 1.1 (or 2.4) for A = 0.2. 



In terms of Wfih„(t) and Jflow (Oi ^he particle number 
in the film, written as -/Vfiini(i), changes in time as 



dt 



iVfilm(i) = W^filmW + Jflow (i). 



(4.13) 



Using the edge velocity fg = dre/dt, we also obtain 



-iT^fiim(i) = 2TTrerene£f, 



(4.14) 



since the film thickness is fixed in our case. In Fig. 10, 
we plot Wtotit), Wmmit), and Jflow(i) vs t for A = 0.002 
and 0.2. In an early stage (t < 1.5 x 10^ for A = 0.002 
and i < 2.6 X 10^ for A = 0.2), Wtot(i) is larger than 
Wfiimit), where condensation occurs on all the interfaces. 
Afterwards, the reverse relation Wtotit) < Wflim(i) holds, 
where evaporation occurs in the droplet body r < rth. 
We also notice Wmmit) > Jflow (*) for t > 1000 for these 
two A. This means that the film extends mainly due to 
condensation near the film edge except in the early stage. 
For example, at t = 1000 (or 4000), we have the edge 
velocity fg = 0.04 (or 0.012) and the gas velocity Vg = 
0.012 (or 0.0064) near the edge in the case A = 0.2. The 
fluid velocity vg in the film is 0.015 (or 0.005) at r == rth. 
These values surely yield Wfn^it) ^ Jflow (0 at t = 1000 
and Wflini(t) ~ 3Jfiow(i) at f = 4000 in accord with their 
curves in the right panel of Fig. 10 and are consistent with 
Eqs.(4.13) and (4.14). Thus, condensation near the film 
edge can be the dominant mechanism of the precursor 
film growth, as originally expected by Hardy [l|, Q . 

We next estimate the gas velocity Vg near the edge. 
The heat flux is of order XiiTp — Tu,)/-^/ there, where 
Tp is the peak temperature and Af is the liquid thermal 
conductivity. It balances with the convective latent heat 
flux ^ UgToAsVg in the gas, where Ug is the gas density 
and As is the entropy difference per particle. Therefore, 



MTp ~ T^)/i£fngToAs) 
(Tp - T^:)neVQ/iTong£f), 



(4.15) 



where we set A^ = kBi^ofie and As(= 2.1kB here) in the 
second line. For example, in the upper plate of Fig. 4 at 
t = 1000 we have Vg = 0.014, while the second line of 
Eq.(4.15) gives 0.012 with If/l - 5 and ni/ug ~ 5. 



V. SPREADING AND EVAPORATION ON A 
HEATED SUBSTRATE 

Next, we present simulation results of a heated liquid 
droplet in the complete wetting condition, where T^ is 
increased above Tq = 0.875Tc at i = with A = 0.2. 
The other parameter values are the same as those in the 
previous section. The preparation method of a droplet is 
unchanged. Then a precursor film develops in an early 
stage (at least for small Ti^ — Tq), because of the complete 
wetting condition at <I>i = 0.061 (see Eq.(3.4)). A new 
aspect is that evaporation dominates over condensation 
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FIG. 12: (Color online) Mass flux J{r, t) through the interface 
in Eq.(4.6) for the weakly heated case T^ = O.SSSSTc corre- 
sponding to (a) in Fig. 11. Evaporation takes place in the 
region J < except close to the edge. Here J decreases in 
time. A steady two-phase state is approached at long times, 
where J is nonvanishing only near the edge. Inset: tempera- 
ture in color and velocity represented by arrows at f = 20000. 



FIG. 11: Time evolutions of the edge position re(i) divided by 
L (left) and the particle number in the droplet Nt(t) divided 
by Ni{q) (right) for A = 0.2. The temperature T„ at the 
solid bottom was raised at i = from O.STSTc to (a) O.SSSSTc 
(top), (b) O.SSSTc (middle), and (c) CSOOTc (bottom). The 
fluid tends to a steady two-phase state in (a), where the liquid 
region assumes a pancake thin film. The liquid evaporates to 
vanish on a time scale of lO'' in (b) and 10* in (c). 



with increasing T^ — To > 0. The experiment by Guena 
et al.^^ corresponds to this situation (see Section 1). 

In Fig. 11, we show the edge position re(i) and the par- 
ticle number in the liquid Ne[t) as functions of t for three 
cases (a) T„ = 0.8855Tc, (b) 0.888Tc, and (c) CSgOT^. In 
the weakest heating case (a) with T^ — Tq = O.OlOSTc, 
re{t) and Ni(t) tend to constants at long times, where a 
thin pancake- like film is realized with radius ~ 0.5L and 
thickness ~ 4£ in a steady state. For higher T^^ , evapora- 
tion dominates over condensation and the liquid region 
eventually disappears. Thus, if T^ — Tq exceeds a critical 
value, a liquid droplet has a finite lifetime due to evap- 
oration even in the complete wetting condition. From 
Fig.ll, this lifetime is of order 10^ at T^ - To = O.OlSTc 
in (b) and 10'* at T^ - Tq = 0.020Tc in (c). 

In Fig. 12, we show the mass flux through the interface 
J{r,t) defined in Eq.(4.6) in the weakly heated case (a) 
Tu, = 0.8855Tc. Its negativity implies evaporation. In 
the region far from the edge, evaporation is marked in 
transient states (i < 4000), but it tends to vanish at long 
times. We can also see the region of positive J with width 
of order 10 near the edge (re — 10 < r < re), where the 
film is still flat and the angle 6 in Eq.(4.7) is nearly tt/2. 
In Fig. 12, however, we do not show J just at the edge 
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FIG. 13: (Color online) Mass flux through the interface J(r, t) 
in Eq.(4.6) for the highest heating case Tw = 0.890Tc corre- 
sponding to (c) in Fig.ll. Here evaporation takes place over 
the whole surface and the liquid disappears at £ = 11000. The 
inset displays the fim shape at i = 3000, where Ne{t) is half 
of the initial value. 



(r = re{t) and < z < £f), where 6 changes from 7r/2 to 
zero in the z direction and evaporation occurs (J < 0). 
As a balance of condensation and evaporation in these 
two regions, the total condensation rate Wtot in Eq.(4.7) 
tends to vanish at long times, while there is no velocity 
field in the region r < r^ ~ 10. In the inset of Fig. 12, the 
velocity field around the edge is displayed at t = 20000, 
where the maximum gas velocity is Vg = 1.1 x 10^'^£/tq. 
In Fig. 13, we show J{r, t) at several times in the high- 
est heating case (c) T^ = O.SOOT^. In the whole surface, 
J is negative and evaporation ocuurs. For t < 40000 
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FIG. 14: (Color online) Top; Temperature T in the region < 
r/L < 0.4 and -0.1 < z/H < 0.7 at t = 1000 in the highest 
heating T^ = 0.890Tc corresponding to (c) in Fig. 11. Bottom: 
temperature in color and velocity by arrows at the same time 
in the same run. A liquid film supports a large temperature 
gradient and evaporation occurs all over the surface, while the 
temperature above the film is nearly flat due to a gas flow. 



evaporation is strongest at the fini center. At long times 
{t — 8000 here), it becomes weakest at the film center. 
Figure 14 is produced by the same run. It gives a bird 
view of the temperature and a snapshot of the velocity 
field in the vicinity of the film edge at i = 1000. We can 
see a steep temperature gradient within the film, which is 
much larger than in the gas, leading to a strong heat flux 
from the solid to the film. In this manner, evaporation is 
induced all over the surface and is strongest at the film 
center in the early stage. It is remarkable that the tem- 
perature gradient nearly vanishes in the gas region above 
the film away from the edge, where heat is transported 
by a gas flow. We can also see a significant temperature 
inhomogeneity in the solid part in contact with the film. 



VI. SUMMARY AND REMARKS 



boundary conditions. The condensation rate on the inter- 
face is a result and not a prerequisite of the calculation. 
We have also assumed that the substrate wall has a finite 
thickness Hy^ and the solid temperature obeys the ther- 
mal diffusion equation, whereas an isothermal substrate 
is usually assumed in the literature. The temperature T^, 
at the solid bottom z — —H^ is a new control parameter 
in our simulation. Cooling (Heating) the fiuid is realized 
by setting T^j lower (higher) than the initial fluid tem- 
perature To- We give salient results in our simulation. 

(i) In the cooled and non-cooled cases with T^, < Tq, 
a precursor film with a constant thickness has appeared 
ahead of the droplet body. Here the liquid volume has 
increased in time due to condensation on the film sur- 
face. In an very early stage, the piston effect comes into 
play due to sound propagation [53, [SJ]- At long times, 
the condensation rate has become localized near the film 
edge and the film has expanded dominantly due to con- 
densation. As a result, a hot spot has appeared near the 
film edge due to the latent heat released. 

(ii) At a critical value of T^ slightly higher than Tq, we 
have realized a steady-state thin liquid film, where con- 
densation and evaporation are localized and balanced at 
the edge. For higher T^,, evaporation has dominated and 
the liquid region has disappeared eventually. This life- 
time decreases with increasing T^ — Tq. For a thin film, 
evaporation has appeared all over the film surface upon 
heating. In our previous simulation for one-component 
fluids 22| evaporation of a thick droplet was mostly local- 
ized near the contact line in the partial wetting condition. 

We give some critical remarks. (1) If the mesh length 
Aa; = £/2 is a few A, our system length is on the order of 
several ten manometers and the particle number treated 
is of order 10^ (see Sec.IIIC). Our continuum descrip- 
tion should be imprecise on the angstrom scale. Thus 
examination of our results by large-scale molecular dy- 
namics simulations should be informative. We should 
also investigate how our numerical results can be used 
or modifled for much larger droplet sizes. (2) In future 
work, we shoud examine the role of the long-range van 
der Waals interaction in the wetting dynamics. As is 
well-known, it crucially influences the film thickness [l|- 
(3) We should also include the slip effect at the con- 
tact line in our scheme [50, l5i|. (4) We should study 
the two-phase hydrodynamics in fluid mixtures, where a 
Marangoni flow decisively governs the dynamics even at 
small solute concentrations [19, 23, ^1 . 



For one-component fiuids we have examined spreading 
of a small droplet on a smooth substrate in the complete 
wetting condition in the axisymmetric geometry. In the 
dynamic van der Waals theory [4l|, |42|, we have inte- 
grated the entropy equation in Eq.(2.14) together with 
the continuity and momentum equations. This method 
may remove artificial flows around an interface [52|. In 
our phase field scheme, we need not introduce any surface 
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